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Abstract 

We have computed by a Monte Carlo method the fourth virial coefficient of free anyons, 
as a function of the statistics angle 6. It can be fitted by a four term Fourier series, 
in which two coefficients are fixed by the known perturbative results at the boson and 
fermion points. We compute partition functions by means of path integrals, which we 
represent diagrammatically in such a way that the connected diagrams give the cluster 
coefficients. This provides a general proof that all cluster and virial coefficients are finite. 
We give explicit polynomial approximations for all path integral contributions to all cluster 
coefficients, implying that only the second virial coefficient is statistics dependent, as 
is the case for two-dimensional exclusion statistics. The assumption leading to these 
approximations is that the tree diagrams dominate and factorize. 

1 Introduction 

Anyons are identical particles in the plane characterized by a statistics phase angle 9 0-|3|. 
They may be thought of e.g. as bosons with a "statistics interaction" between two particles 
which is an angular vector potential proportional to 6/r at a relative distance r. Since this is 
apparently a long range interaction, it is not immediately obvious that the cluster and virial 
coefficients should be finite in the thermodynamic limit. 

We work out here a path integral representation for the cluster coefficients of anyons, 
which shows that they are finite. From this point of view the statistics interaction has short 
range, thus the important property is not the 1/r dependence of the vector potential but 
rather the pointlike nature of the flux. Its range is temperature dependent, however, because 
it is effective when the particle paths wind around each other, and each path in the path 
integral represents Brownian motion of a particle in the plane, covering an area inversely 
proportional to the temperature. 

The path integral representation for the cluster coefficients is in fact quite general, and can 
be applied to anyons in two dimensions, as well as to bosons and fermions in any dimension, 
interacting by general scalar and vector potentials. The same argument for finiteness holds 
in general, when the interaction range is sufficiently short. 

We have computed numerically the fourth virial coefficient of free anyons, by Monte Carlo 
evaluation of the four-particle path integrals. We find that it is very nearly constant, i.e. 
nearly equal to zero for all values of 6. Thus the anyon system is an approximate realization of 
Haldane's so-called exclusion statistics [Q], characterized by a continuously variable parameter 
g, for which only the second virial coefficient depends on g, in two dimensions The 



correspondence between the two kinds of fractional statistics is given by the relation 

g = l-{l-af, (1) 

where a is the periodic function of 9 defined in eq. (|^) below. This correspondence is only 
approximate, in fact it is known from perturbation theory that the higher virial coefficients 
of anyons all have a second order variation with 6 at the boson and fermion points |^12|. 
Our results are consistent with perturbation theory, within the precision obtained. 

All observable properties of anyons must be periodic functions of 6 with period 27r. Energy 
eigenvalues and eigenfunctions are analytic functions of 6, except that some are non-analytic 
at = 0, varying like \9\ rather than 0. Hence the partition functions and all thermodynamic 
quantities derivable from them will be analytic functions of 0, even at the fermion point 6 = tt, 
but generally not at the boson point 9 = 0. 

In the absence of an external magnetic field, the theory is both time reversal and parity 
invariant if each of these transformations is defined so as to include a change in sign of 9. 
It follows that energy eigenvalues and thermodynamic quantities, as functions of 9, must be 
symmetric about 9 = 0, hence they are functions of the quantity a{9) defined by 



a{9) = ^ 
vr 



for 1^1 < ^ , 



a 



+ 2tt) = a{9) . 
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Note that a is non-analytic in 9 at the boson and fermion points, but any even polynomial in 
a is analytic at the boson point, and any even polynomial in 1 — a is analytic at the fermion 
point. An example is the exact second virial coefficient ||l^, 14, |^, 
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with g given as in eq. (||). Here X = h ^2TT(3/m is the thermal wavelength, depending on the 
mass m and the inverse temperature (3. A polynomial in a which is analytic in 9 both at the 
boson and the fermion point, must be constant, because it is a periodic polynomial in 9. 

The third virial coefficient is analytic at the boson as well as the fermion point, because it 
is "supersymmetric" , i.e. symmetric under the substitution a — > 1 — a ||l^, Being analytic 
everywhere, and periodic in 9 with period tt, it should be expandable as a rapidly converging 
power series in s\v?9. In fact, very precise numerical calculations |17], see also [18-20 
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combination with first and second order perturbative calculations for the complete cluster 
and virial expansions |6-12|, give that 



= X' 



1 sm^9 
36 127r2 



(1.652 ±0.012) 10"^ sin^6' + - 



(4) 



It may be conjectured that all virial coefficients, with the exception of A2, are analytic 
functions of 0. If A4 is analytic, then it must have the form 



Aa = X^ 
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^ ln( \/3 + 2 ) + cos 61 + sin^6' (04 + ^4 cos 9) + 
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where the coefficients of the lowest order terms are fixed by perturbation theory at the boson 
and fermion points. The Monte Carlo results presented here can be fitted to this form with 
two parameters. 



C4 



-0.0053 ± 0.0003 



-0.0048 ± 0.0009 . 



(6) 
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We estimated the per degree of freedom (DOF) to be 3.35 for the two-parameter least 
squares fit, and 30 for the no-parameter fit to the minimal Fourier series with C4 = = . . . = 
0. The error estimation is discussed in Section ^ below. We have assumed that the systematic 
errors can be neglected. 

In addition to our numerical results, we derive the following polynomial approximations 
for the cluster coefficients b^, valid for any A^, 

(7) 

One nice property of these polynomials is that they are analytic functions of 9 at the fermion 
point, as the exact cluster coefficients must be. However, they do not give the correct second 
derivatives at the boson and fermion points, known from perturbation theory, although they 
do give the correct first derivatives. An alternative way to introduce the same polynomials is to 
postulate that the second virial coefficient is given by eq. (^), while all higher virial coefficients 
are independent of a. That is, these are just cluster coefficients for two-dimensional exclusion 
statistics |2^, with the statistics parameter given by (||). The corresponding second order 




diagrams were identified in ref. |11]. 



In Sections ^ and ^ below we will review some basic formulae. In Section ^ we describe the 
Monte Carlo method, emphasizing the difference with respect to ref. [19|, where an external 
harmonic oscillator potential was introduced for regularization. In the method used here, we 
regularize by introducing a finite area A. In Section |5| we give some analytical results, either 
exact or approximate. The Monte Carlo results are presented in Section ^, mostly in the form 
of figures. We conclude in Section ^ with a few comments. 



2 The A/^-anyon partition function 



We consider free particles in two dimensions, except that we need to confine them inside a 
finite area A for purposes of normalization. We might use a square box, but it is simpler to use 
periodic boundary conditions so that there are no reflecting walls. Note that the periodicity 
is only used for normalization, and that when we speak about anyons, in the path integral 
formalism, the only restriction is that the starting points of trajectories should be inside the 
given area. Otherwise they propagate freely in the plane and not on the torus. 
The A^-particle partition function is denoted by Zn{(3). In particular. 
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where we have used an identity from ref. |]2^. Below we will use only the leading term in the 
limit A — > 00, that is, we take 



(9) 



The correction terms for finite A are exponential in A. This formula implies the following 
scaling relation, valid for one free particle in two dimensions. 
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By definition, a partition of is a sequence of non-negative integers, V = {i'i,i'2, ■ ■ •)) 
such that N = X^L^i ^^lL. Let Cn denote the set of all partitions of A^, and let C = U??=o 
and C = U?^=i ^N- In this notation we have that 



oo oo 



vec N=oveCN 1/1=0 1^2=0 ul=o 

A partition V ol N labels a conjugation class in the symmetric group Sn of permutations 
of A^ particles, in such a way that vl is the number of cycles of length L. Thus Cn may be 
identified with the set of conjugation classes in Sn- 

The partition function for A^ identical particles can in general be expanded as a sum over 
partitions of A^, 

Zm{P)= MP)Bvm. (12) 

veCN 

Here B-p is the contribution from the partition V to the free particle bosonic partition function, 

and the effect of any interaction of the particles, including the anyonic "statistics interaction" , 
is described by a correction coefficient F-p. 

We consider here the case of free and non- interacting anyons, with an anyon phase angle 
9. In the path integral representation of the partition function, F-p can be interpreted as the 
generating function for the probability distribution of winding numbers |23[| , 

00 

Fp = Fpip,e)= Pr{P,Q) expi-i9Q) . (14) 

Q=— 00 

P-p{P, Q) is the probability of the total winding number Q, given the partition V and given 
the distribution of paths valid for free bosons at the inverse temperature /3. 

Further on, we will label a partition by the set of its addends [so V = (2, 1, 0, . . .), which 
is 4 = 2 + 1 + 1, will be labelled 211]. By eqs. (p^)-(p^) we have in particular that 

Z2iP) = ^FuZi{pf + ^F2Zi{2P) , 

ZsiP) = I Fill ZiiPf + i F21 Zi(2/3)Zi(/?) + 1^3 Zi(3/3) , (15) 
b 2 6 

Zm = ^ Fun Ziipf + ^ F211 Zi(2/?)Zi(/3)2 + i F22 Zi{2pf 
+ i F31 Zi(3/3)Zi(/3) + i F4 Zim . 
Using the scaling (p^), and writing Zn for Zi\i{(3), we get 
Z2 = - Fii Zx' + -F2Z1 , 

Z3 = l Fill ^1^ + 7 F21 Zi^ + ^ Fs Zi , (16) 
6 4 9 



Za — Fi 1 1 1 Z^ ^ + — F21 1 Zi H F22 Zi H F31 Z^ ^ -\- — 

24 8 32 9 16 

For each partition, Zi is raised to a power which is the number of cycles. 
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3 The cluster expansion 

The grand canonical partition function is 



E{(3,ti) = l+J^z''Zr,{f3), (17) 

N=l 

where z = exp{Pfi) is the fugacity and n is the chemical potential. The pressure P is given 
by the relation 

= ^ = E bNz'' , (18) 

^ N=l 

where 6jv is the iV-th cluster coefficient. An immediate consequence is that 



S = n eMAhNz"") = E n ^ , / ' (19) 
iv=i reCL=i ^' 

and hence, 

2«=2:n^- (20) 

We are more interested in the inverse relation, which follows from the expansion 

lnE = ^L^Ij:z^zA = E(-ir-V -l)!n^;T^- (21) 

u=l \L=1 I PgC L=1 

Here v = z^(P) = X^lLi is the total number of cycles in the partition V ■ This gives the 
cluster coefficients in terms of the A^-particle partition functions, 

Ah^= E (-l)-i(z.-i)!n^. (22) 

PeCjv L=i 

In particular. 





= Zi 




Ab2 


= Z2 


' z,^ 

2 ' 


Ab3 


= Z3 


- Z2Z1 + 


Ab4 


= Z4 
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(23) 



zl 

3 ' 

^ + Z2Z^'-^ 
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In general we may write 



AbN = ZN + ---= E i^pSp + • • • = ^1 E G^T' ^ , (24) 
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where we have introduced new coefficients 

Gv = {Fv + -- ■)Zi''-^ . (25) 

The "• • •" in the last formula represents a sum of terms that are products of "-F" coefficients. 

One main point of introducing the "G" coefficients is that they tend to a finite limit in 
the thermodynamic limit yl — > oo, as we will prove below. G-p is the "connected part" of F-p 
for any partition V. The concept of connectedness will also be made more precise below. 

With eqs. (I)-®, we get that 

X'h^ =Y.Gv\{ , (26) 

and all quantities occurring in this equation are finite in the A ^ oo limit. In particular, 

= ^ + ^ + (27) 

We have that Gi = Fi = 1, Gn = Fn fov N = 2,3,4, .. and 



Gil 


= {Fn - 


l)Zi , 


Gill 


= (Fill - 


- 3Fn + 2)Zi^ , 


G21 


= {F21 - 


F2)Zi , 


Giiii 


= (-^1111 


-4Fiii -3Fii2 + 12Fii - 


G211 


= (-P211 - 


- 2F21 — F2F11 + 2F2)Zi^ , 


G22 


= {F22 — 


F2^)Zi , 


G31 


= {F31 - 


Fs)Zi . 



Q)Zi^ , (28) 



For bosons and fermions the probability generating functions can be factorized as 

00 

Fp=l[ Fl"^ , (29) 



where Fl = 1 for bosons and Fl = (—1)^^^ for fermions. This factorization implies that 
which gives the cluster coefficient 

Thus, for bosons and fermions Gat = F^ = (±1)^"^, while G-p = for every partition V 
containing two or more cycles. 
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4 The Monte Carlo method 



We may use the Monte Carlo method in order to compute numerically the coefficient G-p for 
a given partition V representing a conjugation class in the symmetric group S^v- For that 
purpose we represent G-p as a path integral over all paths inducing one given permutation 
represented by V, 

GrZ^=Mv J P(ri(r),...,rjv(r))exp(^-05P- (32) 

Here T'j(r) is the path of particle j, as a function of the imaginary time r, and S is the free 
particle action in imaginary time, 



N .1)3 j„ / \ 2 



5 = E/ dr 







2 



drj(r) 



dr 



(33) 



We include the Gaussian factor exp(— 5/7i) as part of the integration measure, so that it 
is the integrand g-p alone that represents the interaction of the particles, and we include a 
normalization factor Afp so that Gp = Zi^"^ \i gp = 1 identically. Note that Np is then 
finite (i.e. A independent), since the path integral is proportional to Zi^ when gp> = 1. Note 
also that this path integral representation is actually very general, and can be applied to any 
iV-particle system with (short range) interactions in any dimension, not just to the iV-anyon 
system considered here. 

To see what the integrand g^p looks like in our case, let us take the partition 2+1+1 of 4 as 
an example. A closed path in the four-particle configuration space interchanges the positions 
of two particles, say particles 1 and 2, and takes the remaining two particles back to their 
starting points. The total winding number Q is the sum of six pairwise winding numbers, 

Q = Qi2 + {Qi3 + Q23) + {Qu + Q24) + Q34 . (34) 

Note that Q12 is an odd integer and (534 an even integer (remember that the winding numbers 
are defined such that a complete revolution corresponds to the winding number 2), whereas 
<3i3) Q23, Qu, Q24 are in general non-integer, because particles 1 and 2 do not return to their 
starting positions. However, the sums (5(i2)3 = Q13 + Q23 and (5(i2)4 = Qi4 + Q24 are even 
integers. Hence Q is an odd integer. Let / be any subscript, and introduce the notation 

ei = l + fi = exp(-i^Qj) . (35) 

In order to compute the coefficient G211Z1 = (F211 — 2F21 — -F2F11 + 2F2)Zi^ we take the 
integrand to be 

3211 = ei2 6(12)3 6(12)4 ^34 — 612 6(12)3 ~ ^12 6(12)4 " 612 634 + 2ei2 

= 612 (/(12)3 /(12)4 /34 + /(12)3 /(12)4 + /(12)3 /34 + /(12)4 /34) • (36) 

For example, we compute F2iiZi^ by integrating 

exp{-i9Q) = ei2 6(12)3 e(i2)4 634 , (37) 
and we compute 2^21^1'^ by integrating 

exp(-i6l(Qi2 + <3(i2)3)) + exp(-i6'((5i2 + Q(i2)4)) = 612 6(12)3 + ei2 6(12)4 • (38) 
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The second line of eq. (36) may be represented diagrammatically as 



^211^1 




+ 




+ 2 



(39) 



The particles are represented as points (filled circles). The two-cycle is represented by eu in 
the integrand and by a circle connecting two particles in the corresponding diagram. Each 
factor // in the integrand is drawn as a single straight line in the diagram. Note that we 
should draw labelled graphs to represent the four terms in eq. (|36|). But since the value of 
a graph is independent of the labelling, it is more natural to draw unlabelled graphs and 
include instead integer coefficients counting the number of ways each graph can be labelled. 
Hence the factor 2 in front of the last graph. 

In a similar way we find the diagrammatic representation 



GiiiiZi 



1=^+6 +12 ]/l 



+ 3 



+ 4 



+ 12 



(40) 



The coefficient in front of each diagram is again the number of inequivalent ways of labelling 
the nodes of the graph. We may also write 



G22Z1 



G31ZI 



(41) 



We see that only connected diagrams contribute to the cluster coefficients. It follows that 
the latter are finite in the limit j4 — > 00. Indeed, any path gives a non-zero contribution 
to the path integral represented by some diagram only if for every line in the diagram, the 
corresponding winding number is non-zero. The probability for this to happen for a connected 
diagram goes to zero as {X'^/AY~^ when A — > 00, since every L-cycle path gives a Gaussian 
distribution of points which essentially covers only a finite area, proportional to A'^. Here v 
is the number of cycles, and — 1 is the minimum number of links in a connected graph with 
u nodes. The factor A~^^^ cancels exactly the divergence of the factor Zi^~'^ included in the 
definition of G-p, eq. (p8|). 

The general meaning of the relations between the F and G coefficients should now be 
obvious. F-p is a sum of both connected and disconnected diagrams, whereas G-p is the part 
of the sum including only the connected diagrams. For example, the relation 



F211Z1 — G211Z1 + 2G21G1Z1 + G2G11Z1 + G2G1G1Z1" 



(42) 



which follows from (|28D, is represented as 



-p2ii^i" 




+ 




G211Z1 



+ 2 



+ 



o 



+ 



o 



(43) 



2G2iGiZi'^ G2G\\Z}^ G2G1G1Z1 



It is the last term that dominates in the thermodynamic limit, but it is G211 only that 
contributes to the cluster coefficient. Thus, as usual, the grand partition function is a sum of 
all diagrams but the thermodynamic potential is a sum of connected diagrams 



The Monte Carlo method consists in generating random paths according to the Gaussian 



distribution of paths valid for bosons |19|. Each four-particle path is closed over the imaginary 
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time interval hp, in the sense that the final configuration is identical to the initial one, but 
with the particle positions interchanged by a permutation belonging to the class V C Sat. 
Consider the partition 2 + 1 + 1 = 4, as in the example above. Then particles 1 and 2 should 
interchange positions, while particles 3 and 4 should return to their starting points. We take, 
arbitrarily, the starting point for the path of particle 1 to be at the origin, this is then also 
the ending point for particle 2. Equivalently, it is the ending point for particle 1 over the 
imaginary time interval 2hp. The starting point for particle 2, equal to the position of particle 
1 after half the imaginary time interval 2hp, can then be generated according to a Gaussian 
distribution around the origin. The starting and ending point for particle 3 is generated 
according to a flat distribution inside a square area A centered on the origin. Similarly for 
particle 4. 

For each four-particle path generated we count the winding numbers Qu, Q(i2)3) Q(i2)4) 
Q34 and increment a histogram n{Q) in the following way. We compute the total winding 
number Q and add 1 to n{Q), this takes care of the integrand eu 6(12)3 6(12)4 634- We subtract 
1 from n{Qi2 + Q(i2)3)> ™ order to take care of the integrand —6126(12)3- Similarly, we 
subtract 1 from n{Qi2 + Q(i2)4) and from n{Qi2 + Qsa), and we add 2 to n((5i2)- Finally, 
G211 is the Fourier transform of the histogram n{Q), multiplied by the normalization factor 
Zi^/n, where n is the total number of four-particle paths generated. The net contribution to 
the histogram vanishes if more than one of the three winding numbers (5(12)3, Q{i2)4 &iid Q34 
is zero, and this is what ensures a finite limit as ^ — > 00 for the computed G211. 

5 Exact and approximate polynomials 

The first cluster coefficient, with our definition, is 61 = The exact result for the second 

cluster coefficient of free anyons is 

X% = ^ + ^ = lil-af-l. (44) 

Since Gu is even and G2 is odd under the substitution a ^ 1 — a, this implies that 

Gil = a(a - 1) , G2 = F2 = 1 - 2a . (45) 
It is further known that the third virial coefficient, 

A3 = -2-^ +4^, (46) 
bi &i 

is even under a — > 1 — a. The odd part of —263/61'^, which is — A^G2i/2, therefore has to 
cancel the odd part of 462'^/6i^, which is A^GiiG2. This condition gives another exact result, 

G21 = 2F2 Gil . (47) 



One further result |25|, which is exact according to the perturbative calculation of ref. |2^, is 

L-l 



n(i-^)- m 



k=l 



Our Monte Carlo results for single cycles of length L < 4 are consistent with this formula, 
which is a check of our Monte Carlo calculation as well as of the perturbative calculation. 
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Another way to compute G2 = F2, Gn and G21 is to impose an external harmonic 
oscillator potential to make the energy spectrum discrete, and then take the zero frequency 
limit. Only energy levels depending linearly on a contribute to these three quantities. The 
computation of Gl = Fl for L > 2 is much more non-trivial, and in fact the only known 
method is perturbation theory, because also states with non-linear a dependence contribute. 

This is about as far as one can get with exact results. However, in the diagrammatic 
expansions shown above, one may argue quite generally that the tree graphs are expected to 
dominate, because every additional line in a diagram represents another factor of the type 
// = exp{—i9Qj) — 1 in the integrand, with Qj an even integer. This factor vanishes when 
Qi = 0, which will happen with a certain probability which is definitely non-zero, and even 
if it does not vanish it will often have an absolute value smaller than 1. Furthermore, one 
may argue that the path integral represented by a tree graph should approximately factorize 
in the same way as its integrand. These two assumptions, of tree diagram dominance and 
factorization, lead in a not entirely trivial way to the following polynomial approximation for 
the general coefficient G-p, 

00 

Gv^Gv = N-'-^ Gu"-^ n (^^i)'^ • (49) 

L=l 

There is a factor for every cycle of length L, a factor L1L2G11 for every single line 
connecting two different cycles of lengths Li and L2 (each L-factor counts the number of 
ways the line can be connected to the cycle), and there is a sum over all possible ways 
to connect the cycles into a tree graph. It is perhaps not obvious how this leads to eq. (^9|); 
a simple way to understand the connection is by looking at low order examples: Consider 
the case of 3 cycles of lengths Li, L2, and L3. They can be connected to a tree graph in 3 
possible ways. This gives a coefficient 

FliFl^Fl^ (L1L2G11 L2L'iGii + L2L3G11 L3L1G11 + L3L1G11 L1L2G11) 

3 

= (Li + L2 + ^3) Gn^ WUFl^, 

i=l 

which agrees with eq. (H) since Li + L2 + L3 = A^. We should point out that eq. (|^) was 
first derived empirically as an approximate representation of the Monte Carlo results. 

Special cases where these polynomial formulae are exact, as already mentioned, are the 
cyclic coefficients Gl = Fl, as well as Gn and G21 = 2F2G11. In the three-particle case 
there is one approximate polynomial, 

Gill = 3Gii2 . (50) 

The four-particle approximate polynomials are: 

Giiii = IGG'ii'^ , G211 = 8F2Gii^ , G22 = 4F2^Gii , G31 = 3F3G11 . (51) 

The polynomial approximations for the G coefficients imply the polynomial approxima- 
tions for the cluster coefficients given in eq. (0), and those imply in turn that the virial 
coefficients are independent of the statistics, except for the second coefficient (see Appendix). 
These polynomials (0) are characteristic of two-dimensional exclusion statistics with statistics 
parameter (|l]). 
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6 Monte Carlo results 



In the Monte Carlo integration one has to choose a suitable value for the dimensionless 
parameter ^/A^, making a compromise between systematic and statistical errors. We want 
the limit Aj)^ oo, but if we take A/X'^ too large, we get no statistics. As we have seen, 
one randomly generated four-particle path contributes to the numerical path integral only if 
the paths of all four particles wind around each other as one cluster. This necessary condition 
implies that the systematic errors for large but finite A/X^ are exponentially small, since the 
maximum distance that one particle wanders away from its starting point, has a Gaussian 
distribution. The exponential convergence for the G coefficients is similar to, although not 
directly related to, the convergence for Zi. 

Exponential convergence means that one should use one single value of A/X"^, as large as 
practically useful, rather than use two or three values and extrapolate. We used the following 
values, 

partition 1+1+1+1: tt^tt^ = 3.183 , (52) 

ZTT U.U5 

partition 2+1+1: ^ I nn = 5-305 , (53) 

ZTT U.Uo 

partitions 2+2 and 3+1: = 15.915 . (54) 

ZTT 0.01 

Typically, two to five out of 1000 generated four-particle paths gave non-zero contributions 
for these partitions. We have generated 1.7 10^ paths for the partition 1 + 1 + 1 + 1, 16.4 10^ 
for 2 + 1 + 1, 1.5 10^ for 2 + 2, 17.6 10^ for 3 + 1, and finally 397000 paths for the partition 4. 

The computed G coefficients are plotted in the Figures || to ^ as functions of a. In each 
case the polynomial approximation, as discussed in the previous Section, is subtracted, and 
the resulting curve is marked "Re(MC) — polynomial". Because of the statistical errors, the 
Monte Carlo generated curve has also a non-zero imaginary part, marked "Im(MC)", which 
is useful because it indicates the statistical errors in the real part. Since the real part is 
even about a = 1/2 and the imaginary part is odd, or vice versa, depending on whether the 
partition is even or odd, we always plot only the interval < a < 1/2. 

Fig. H shows the computed Gnu with the polynomial 166*11^ = 16(a(a — 1))'^ subtracted. 
The curve marked "fit" is mostly empirical, and is given by 

3 

fit = 5- a(a — 1) sin^(a7r) . (55) 

The figure shows that this is a perfect fit to the Monte Carlo curve, within the statistical 
uncertainty as indicated by the imaginary part. 

Fig. H shows the computed G211 with the polynomial 8F2Gii^ = 8(1 — 2a){a{a — 1))^ 
subtracted. The curve marked "fit" is partly empirical, but with a coefficient which is chosen 
so as to produce the correct second order derivative at a = p^]. The formula is: 

fit = ^ (1 - 2a) sm\aTT) . (56) 

37r^ 

Fig. I shows the computed G22 with 4^2^011 = 4(1 - 2afa{a - 1) subtracted. The "fit" 
here is 

2 

fit = ln(\/3 + 2) sin ^tt) cos V^t) • (57) 

V37r^ 
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Fig. ^ shows the computed G31 with 3F3G11 = 3(1 — 3a) (1 — (3/2)a)a(a — 1) subtracted. 
Here we have chosen 

/3 

fit = ^ ln(\/3 + 2) sin^vr) cos V^r) • (58) 

Fig. ^ shows the computed G4 = F4 with the polynomial (1 — 4a)(l — 2a)(l — (4/3)a) 
subtracted. The figure supports the claim that the polynomial is exact. 

Fig. ^ shows the computed cluster coefficient, X^b4 with the polynomial A^64 of eq. (Q) 
subtracted. The parabolas given by the second order perturbation theory at a = and a = 1 
are shown. 

Fig. shows the computed virial coefficient, A^/X^. The parabolas given by the second 
order perturbation theory at a = and a = 1 are shown. Also plotted are two Fourier 
series, as given in eq. (^. The curve marked "Fourier 1" has only the two terms required by 
perturbation theory, i.e. C4 = (^4 = . . . = 0, whereas the curve marked "Fourier 2" is a least 
squares fit with the coefficients C4 = —0.0053 and ^4 = —0.0048. 

The error estimation for the least squares fit requires some comment. We neglected the 
systematic errors, assuming that they are small. Since the imaginary part of A4 is zero, 
by definition, we may estimate the P^r degree of freedom (DOF) for any fit to the real 
part of the Monte Carlo data by scaling so that x^/DOF = 1 for a comparable fit to the 
imaginary part. For the no-parameter fit to A4 with C4 = ^4 = . . . = 0, we get x^/DOF = 30 
by comparison with the no-parameter fit of the exact value zero to the imaginary part. The 
same scale factor would give x^/DOF = 1.18 for the two-parameter fit to A4 with C4 and as 
free parameters. However, the two-parameter fit for the real part should rather be compared 
to a two-parameter fit of the same two terms to the imaginary part, and when we do so, 
we get instead x^/DOF = 3.35. Comparing the no-parameter and two-parameter fits to the 
imaginary part, we find a third scale factor such that these two degrees of freedom contribute 
two to x^- The errors of the two fitted parameters C4 and for the real part are defined as 
the changes corresponding to an increase in x^ of one, when this third scale factor is included. 
This gives C4 = -0.0053 ± 0.0003 and ^4 = -0.0048 ± 0.0009, with uncorrelated errors. We 
arrive at the clear conclusion that the minimal Fourier series, i.e. with C4 = ^4 = . . . = 0, is 
excluded by our Monte Carlo results. 

We show Figures |^, ^ and |10| in order to justify partly the polynomial approximations to 
the tree graphs. The quantity plotted is the contribution G^j^l^^ to a tree graph of one single 
line, joining two cycles of lengths (^1,^2), in the approximation that the integral represented 
by the graph factorizes. We claim that this contribution is L1L2G11, to lowest order. Fig. ^ 
has Li = 2 and L2 = 1, and the "fit" which is subtracted from the Monte Carlo data is 

fit = 2a{a - 1) + 1.07(a(a - 1))^ - 0.3(a(a - 1))^ . (59) 

The important information is the coefficient L1L2 = 2 of the term Gn = a{a — 1). Fig. ^ has 
Li = L2 = 2, and the "fit" 

fit = 4a(a - 1) + 5.77(a(a - 1))^ + 2.2(a(a - 1))^ (60) 

has been subtracted. Fig. |lO| has Li = 3, L2 = Ij and we have subtracted 

fit = 3a(a-l) + 2.46(a(a-l))2 . (61) 
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The (unphysical) sharp peaks near a = in these figures are due to a deficit of high winding 
numbers in the Monte Carlo simulations. 

The last two figures show the splitting of Gm into a tree graph and a triangle graph, 
as the simplest possible example of the contributions of separate graphs. We Monte Carlo 
generated 40.710^ three-particle paths at A/X'^ = l/(27r0.03) = 5.305. Fig. 11 shows the tree 
graph contribution G^ff, after the following "fit" has been subtracted 

fit = 3(a(a - l)f - (l/(47r2)) sm^an) - 1.41(a(a - 1))^ + 1.8(a(a - 1))^ . (62) 



Fig. 12 shows the triangle contribution G*"^^"^^'^, after subtraction of 

fit = 1.41(a(Q - 1))^ - 1.8(a(a - 1))^ . (63) 

This example should give some general idea of what errors we make by the two approximations 
of neglecting loop diagrams and representing the tree graphs by the polynomials of eq. (^9|) . 
We find no apparent simple mathematical form for the correction terms. For the tree graph 
alone the leading correction to the polynomial 3(a(a — 1))^ is of order at the boson point, 
and (a — 1)^ at the fermion point. The sin^ form of this correction, as given in eq. (|62|), 
reproduces the well known third virial coefficient. In general, we would expect corrections 
of order a" and (a — 1)" from graphs with n links, exactly as seen in the present example. 
However, if we consider not the perturbative orders in a and a — 1, but rather the magnitudes 
of the maximal corrections from the tree and triangle graphs, we find that they are comparable. 

These considerations suggest that the splitting into graphs does not simplify the problem 
of calculating the cluster coefficients. One would like for example to develop a perturbation 
theory for calculating the contributions from separate graphs, but at present we have no such 
theory. 



7 Conclusion 

The main results of this paper are the following. The fourth virial coefficient of anyons as 
computed by the Monte Carlo method is rather well described by a Fourier series consistent 
with second-order perturbation theory and deviating very little from zero. Our diagram 
analysis shows in general that all cluster and virial coefficients are finite, and gives, in a 
certain approximation, a direct correspondence with exclusion statistics. In particular, in 
this approximation all the virial coefficients starting from the third are constant. 

The above approximation, as well as the fact that also the exact third and fourth virial 
coefficients apparently exhibit quite a regular behavior with a, seem to hint that the inter- 
polation between Bose and Fermi statistics is rather smooth. This is in strong contrast to 
e.g. the anyon superconductivity phenomenon (non-analytic in 9) believed to take place at 



low temperatures |28-31|. However, since the dimensionless expansion parameter of the virial 
series is x = A^p, which goes to infinity at zero temperature, it is difficult to predict the low 
temperature behavior from the first few terms of a virial expansion. 

What we can test, on the other hand, is the average field approximation on which 
many of the low temperature predictions are based. In this approximation one replaces the 
anyon statistical flux with an average magnetic field B = vp^Q (with = h/e the elementary 
quantum of fiux), and treats the anyons as a collection of non-interacting bosons {ly = 9/tt) or 
fermions {v = O/ir — l) in the average magnetic field. This approximation should work equally 



13 



well at high and low temperatures. Thus, we may compare the resulting virial expansion with 
our results. Let f±{p,B) be the Helmholtz free energy of non-interacting bosons/fermions 
in a magnetic field B. Then f±{p,vp^o) is the free energy of anyons in the average field 
approximation. This leads to a virial expansion 

^„,2 1 2 l + 3i/2 1/2 1 - I00z^2 ^ 5j,4 

/3PA2 = X =f - x2 + x^±—x^ + ■■■ (64) 

^ 4 36 16 3600 ^ ' 

Not unexpectedly, eq. ( p^ ) does not reproduce fermions (bosons) correctly when expanded 
about bosons (fermions). Implementation of such behavior requires corrections which have to 
be put in by hand. If this is done (by e.g. a minimal Fourier series), it seems that the average 
field virial coefficients have a magnitude, and variation with 6, which up to 4th order in x is 
qualitatively rather similar to our results. Since the average field approximation predicts a 
thermodynamics with very interesting low temperature behavior, our results do not rule out 
that this will happen in the exact model as well. 
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A Polynomials for cluster and virial coefficients 

In this appendix we will prove the result that the polynomial approximation (0) for the cluster 
coefficients is equivalent to a virial expansion which is the same as for the two-dimensional 
free non-relativistic Bose gas, except that the second virial coefficient is given by eq. d). We 
will also prove that the polynomial approximation (^) for G-p implies eq. (0). 

For simplicity we fix the temperature and choose units such that /3 = A = 1. Thus, e.g., 
the fugacity is 2; = e^. We make use of the expansions 

dp p dp 



We also define 



00 N-l / ^^^^ 00 ^ r\N-l 



Ar=l k=l ^ ' N=\ 



which is the density corresponding to the cluster coefficients of eq. (0). For the Bose gas, 
with g = 0, we have p = — ln(l — z). Shifting the second virial coefficient by an amount 
AA2 = 5/2 then gives 

H = ln{l- e'f) + gp . (67) 
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For every g > and every /x, or for g = and every /U < 0, this equation clearly has a unique 
solution p > 0. We want to prove that the solution is p = PgifJ-)- 
For this purpose we rewrite eq. ([67|) as 



oo 



p = -ln(l-.e-^0 = E - 



-ngp 



(68) 



n=l 



and apply the following theorem due to Lagrange |27]: The equation p = f{p) has the solution 



r=0 



This gives 



oo -j^ oo 

/'= E E 

M=l ni=l 



E 



^.niH hriM / d \ ^^^1 

- I I 



-(niH VnM)gr 



5: (-A^<7)''-^ 



r=0 



AT 



Ar=l M=l 



where 



Cat, 



1 



M 



oo r 

Oni+---+nM,N 



E--- E 

ni=l njv/=l 



What we need to show is that 



N 



M=l 



(_l)iv-i fNg 



It is straightforward to show that 

oo N 



-gln(l-z) 



N=l M=l 



and hence, 



N 



E 9^'-' Cn, 



Ng \ N 



i=E( 

N=l 



\N 



N (-9 
N 



M 



M=l 



9 \N 



Substituting g — > —Ng we get eq. (|72|), completing the proof. 
We next turn to the cluster coefficients 



(69) 



(70) 



(71) 



(72) 



(73) 



(74) 



V^Cn 



(75) 



given by the polynomial approximation in eq. (^9|). We want to prove that b' 
We may rewrite the above formula as 



AT — bN- 



b'N=Y: 



■ E ■ ■ ■ E ^nr+-+n,,N n 
ni=l n^=l 



• 1 "^j 



(76) 
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To evaluate p = Y^^=i ^^'n insert eq. ([76|), interchange the summation order of N and 

u and use the relations N'^-^z^ = {d/dpY-^z^ and Y.n=i ^"^n/^ = pa{^J)■ We find 



u=l 



(77) 



r=0 



By the Lagrange theorem, eq. (|7^ ) is the solution to the equation p = pa{fJ. + Gup), which, 
as we saw above, is equivalent to 



p + Giip = In (l — e ^) + ap . 



(78) 



This is precisely eq. ( |67| ) with g = a — Gn = 1 — (1 — a)^, which means that 6^ = bj\[ with 
^AT as given in eq. (0). 
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Figure 1: Gim — 16(a(a— 1))^ as a function of a. The imaginary part is plotted to indicate the 
statistical uncertainty of the real part of the Monte Carlo data. Only the interval < a < 1/2 
is plotted, because of the (anti) symmetry about a = 1/2 The curve marked "fit" is given in 
eq. dH). 
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Figure 2: G211 — 8(1 — 2a){a{a — 1))^ versus a. The curve marked "fit" is given in eq. (|56|). 
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Figure 3: G22 — 4(1 — 2a) a{a — 1) versus a. The "fit" here is given in eq. (| 
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23 




Figure 6: The fourth cluster coefficient minus the polynomial of eq. A^(64 — 64), as a 
function of a. Also shown are the parabolas given by the second order perturbation theory 
at a = and a = 1. 
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Figure 7: The fourth virial coefficient, A^^/X^, as a function of a. Also plotted are the 
parabolas given by the second order perturbation theory at a = and a = 1 , and two different 
Fourier series, as given in eq. (|5|). The curve marked "Fourier 1" has C4 = (i4 = . . . = 0, 
whereas "Fourier 2" is the least squares fit with C4 = —0.0053, = —0.0048. 
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Figure 8: G^Ji^ - 2q:(q! - 1) - 1.07(a(a - l))^ + 0.3(a(a - 1))^ versus a. 
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Figure 9: G^J^ - 4q:(q! - 1) - 5.77(a(a - l))^ - 2.2(a(a - 1))=^ versus a. 
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Figure 10: G|f^ - 3a{a - 1) - 2.46(a(Q: - l))^ versus a. 
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Figure 11: The tree graph contribution G^'ff minus the "fit" given in eq. (|62|). 
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Figure 12: The triangle graph contribution G^l^'^^^ minus the "fit" given in eq. (| 
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